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A new implementation of the GW approximation (GWA) based on the all-electron Projector- 
Augmented- Wave method (PAW) is presented, where the screened Coulomb interaction is computed 
within the Random Phase Approximation (RPA) instead of the plasmon-pole model. Two different 
ways of computing the self-energy are reported. The method is used successfully to determine the 
quasiparticle energies of six semiconducting or insulating materials: Si, SiC, AlAs, InAs, NaH and 
KH. To illustrate the novelty of the method the real and imaginary part of the frequency-dependent 
self-energy together with the spectral function of silicon are computed. Finally, the GWA results are 
compared with other calculations, highlighting that all-electron GWA results can differ markedly 
from those based on pseudopotential approaches. 
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I. INTRODUCTION 

For many weakly correlated materials, the density- 
functional theoryi(DFT) in the local-density approxima- 
tion (LDA) provides a good description of their ground- 
state properties. However, DFT is not able to describe 
correctly their excited states. Thus, for example, the 
band gaps in the LDA are typically much smaller than 
the experimental values. Quasiparticle (QP) electronic- 
structure calculations beyond the DFT are therefore 
highly desirable. 

The GW approximation (GWA) of Hedin 2 - 3 -, which 
produces a good approximation for the electron's self- 
energy E, enables us to make first-principle QP calcula- 
tions for realistic materials. Thus, the GWA has been 
successfully applied to the calculation of QP electronic 
structures of many kinds of materials^*^. In particu- 
lar, recent success has been achieved on predicting the 
metal-insulator transition of bcc hydrogen^ electronic 
excitations of yttrium trihydride^ as well as the QP elec- 
tronic structure of copper iiS Unfortunately, most of the 
GWA implementations are based on the pseudopoten- 
tial type of approaches together with plasmon-pole (PIP) 
modelsiii*i2iiiii^iii*i& The weakness of these types of cal- 
culations is that the imaginary part of the self-energy is 
not accessible, making it impossible to determine spec- 
tral functions and hence to interpret photoemission ex- 
periments. In addition the PIP approximation is ex- 
pected not to hold for systems with localized electrons. 
Moreover, it has been noticed recentljii^i^ that GWA 
implementations based on pseudopotential methods lead 
to larger and more k-dependent shifts than calculations 
based on all-electron DFT methods, bringing into ques- 



tion the validity of the former approaches. 

However, some attempts have been made to go beyond 
the plasmon-pole approximationAiSiiLi&iiiSSiSl In par- 
ticular, Aryasetiawan has approximately determined the 
screening within the RPA using a linear muffin-tin or- 
bital method within the atomic sphere approximation^ ( 
LMTO-ASA). This method, although fast, approximates 
the space by atomic centered overlapping spheres, thus 
completely neglecting the interstitial region, and hence 
making the reliability of the GW method uncertain. 
Kotani and van Schilfgaarde based their full-potential 
LMTO GW calculation 17 on the work of Aryasetiawan by 
taking into account correctly the interstitial region. Nev- 
ertheless their method is not quite accurate since in their 
implementation they didn't take into account the multi- 
plicity of the same angular momenta for a given principal 
quantum number in the basis set (like simultaneously us- 
ing the 3d and 4d states). Finally, Ku and Eguiluz pro- 
duced selfconsistent and non-selfconsistent QP band gaps 
based on an approximate Luttinger-Ward functional, 18 
the non-selfconsistent results are much smaller than all 
existing GW calculations. Since these results are based 
on a different scheme we have chosen not to discuss their 
method further. On the other hand, several pseudopoten- 
tial have produced GW results without resorting to the 
plasmon-pole approximation. These methods, although 
interesting, use pseudowave functions and hence can only 
determine pseudo-matrix elements of operators, making 
them difficult to justify as quantitative and reliable meth- 
ods for computing QP properties. 

The major purpose of this paper is then to present 
a new implementation of the GWA method using 
the all-electron full-potential Projector Augmented- Wave 
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method (PAW) complete basis set, and without using any 
PIP model for the determination of the dielectric func- 
tion. The screening of the Coulomb interaction is thus 
described in the Random Phase Approximation (RPA), 
avoiding further approximations. 

The paper is organized as follows. In section |H] we de- 
scribe our implementation of the GW approximation. In 
Section UTTl we present our QP calculations for Si, SiC, 
AlAs, and InAs and also for the alkali hydrides com- 
pounds NaH and KH (to our knowledge this is the first 
GWA study of these alkali hydrides). At the end of this 
section we compare and discuss our results with other 
calculations and experiments. 

II. FORMALISM 
A. The PAW method 

The PAW formalism has been well-described 
elsewhere ^S^SliSk so we will not discuss it in this 
paper. The PAW method- 2 - is a very powerful all- 
electron method for performing electronic structure 
calculations within the framework of the LDA. It takes 
advantage of the simplicity of pseudopotential methods, 
but describes correctly the nodal behavior in the aug- 
mentation regions. The selfconsistent calculation of the 
electronic structure is performed using the Car-Parincllo 
method over the occupied states. To determine the 
eigenvalues and eigenvectors of all unoccupied states (up 
to 200 eV above the top of the valence states) needed 
for the GW calculations, we have extracted the selfcon- 
sistent full-potential, constructed and diagonalized the 
PAW Hamiltonian for every irreducible k point in the 
Brillouin zone. 



B. The GW approximation 

1. Quasiparticle energies 

In general, the QP energies E n (k) and wave function 
"0kn(r) are determined from the solution of the QP 
equation 

(T + V ext + WknW + /d 3 r'S(r,r',S„(k))^k„(r') 
= £„(k)V> k „(r) (1) 

where T is the free-electron kinetic energy operator, V ex t 
the external potential due to the ion cores, V% the average 
electrostatic (Hartree) potential, and £ the electron self- 
energy operator. The major difficulty connected with 
Eq. JQ) is finding an adequate approximation for the 
self-energy operator S(r, r', E n (k)). Nonetheless, it was 
shown by Hedin 2 ^ that writing the self-energy as a prod- 
uct of the Green's function and the screened Coulomb 



interaction W yields the successful GW approximation 
for E. In this approximation, both the non-locality and 
the dynamical correlations are included. Assuming that 
the difference X — V xc between the self-energy and the 
Kohn-Sham exchange and correlation potential is small, 
we can use a perturbation theory approach to solve the 
effective QP Hamiltonian H qp 

H qp = Hks + {t- V xc ) (2) 

and determine the QP energies by expanding the real 
part of selfenergy to first order around e n (k) thus making 
the comprison with the PIP models possible 
ReE n (k)) = e„(k) +Z„ k x 

[<tf kn |Re£(r,r',e n (k))|tf kn > - {^k n \V x L c DA (r)\^ kn )} (3) 
where the QP renormalization factor Z nk is given by 

Z„k = [1 - (* kn |^-ReE(r,r',e„(k))|* k „)]- 1 . (4) 

OiO 

This assumption is valid for simple sp bonded materials, 
since it was shown that the QP wave function V'kn and 
Kohn-Sham wave function are almost identical, i.e., 
the QP Hamiltonian H qp is diagonal in the \I/ kn basis 
for simple sp bonded semiconductors^^ We therefore 
assume this behavior for the materials studied in this 
paper. According to this equation, the LDA eigenval- 
ues e„(k) are then corrected by the GW approximation. 
The numerical work is therefore considerably reduced, 
but still computationally demanding. 

In our implementation, we have calculated the Green's 
function only for the valence and conduction states. One 
has then to subtract out only the valence exchange and 
correlation potential in Eq. <|3J). To check the accuracy of 
this procedure, we have also used the so-called Hartree- 
Fock decoupling and have found that the average er- 
ror in the QP energies of Si with respect to the top of the 
valence states is 0.05 eV. The approximation made here is 
the one currently used in all pseudopotential-based GWA 
calculations, making our method compatible with exist- 
ing GW implementations. 



2. Screened Coulomb interaction 

For the calculation of the self-energy, one needs to eval- 
uate the dynamically screened interaction W(r,r',ui), 
which can be rewritten in reciprocal space as: 

W 6 ^ (q)W ) = ^_^g a i G/(q)W )_J_. ( 5 ) 

The symmetrized dielectric matrix ?gg' (q, is defined 
in the random phase approximation (RPA) bySi 
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fi|q + G||q+G 



-]TM G c (k,q)[M G c ,(k,q) 



1 



l 



u + e t ,(k — q) — e c (k) — iS u — e„(k — q) + e c (k) + iS 
I 



(6) 



with the following notation: 

Ma m (k,q) = <vP k _ q „|e-^+ G >|*k m }, (7) 

where v and c denote, respectively, the valence and con- 
duction states, and 5 a positive infinitesimal. The matrix 
elements given by Eq. are evaluated using the PAW 
basis set as described in Ref. 0. 

Most of GWA calculations use a kind of PIP approx- 
imation. This is computationally efficient since one ob- 
tains an analytic expression for the integral in the self- 
energy. It is not clear, however, that this kind of approx- 
imation is valid for describing the QP of different kind 
of materials. It is for this reason that we have chosen to 
avoid the PIP model altogether, to compute the dynami- 
cal dielectric function in the RPA JBJ, and to perform the 

I 



integral of the selfenergy numerically. In our implementa- 
tion, we need to compute ecc (q, w) along the imaginary 
axis and for some real frequencies. This technical point 
will become clearer in the next subsection. 

To reduce the computational cost of the GWA, we use 
symmetry properties. Details about the utilization of the 
symmetry for the static dielectric matrix has been al- 
ready given elsewhere*i*i£*22i2& ; so we just describe briefly 
how to use the symmetry in the case of the dynamical 
dielectric function. For the case of pure imaginary fre- 
quencies, we could safely ignore the broadening factor iS; 
in this case e GG ' (q, iu) is hermitian and we could use the 
symmetry just as in the static case. We can then write 
the symmetrized dielectric matrix as 



e GG <(q,iw) = fcc - 



8tt 



ft|q+G||q + G 



1 



^ E EE Mj£(k,q)[M£ c G ,(k, q )r 



k£BZ„ v,c ReG„ 



- e„(k - q) - e c (k) iu - e„(k - q) + e c (k) / ' 



(8) 



where G q is the little group of the point group G such 
that i?q = q; R being a symmetry operation. The com- 
putational cost is further reduced by noticing that 



e G G'(-Rq,M = ei?-iGfl-iG'(q,^) 



(9) 



for real u>, although the dielectric matrix is not hermi- 
tian, we could use the symmetry by making a decom- 
position into hermitian and anti-hermitian parts of the 
polarizability PGG'(q,w). if we define AGG'tq,^) and 
S G G'(q,w) by 



AcG'(q,w) 



P GG '(q,w) + f GG ,(q,t 



(10) 



and 



then equations JSJ and JJjJ still hold, allowing us 
to perform the same computational tasks as for the 
symmetrized dielectric matrix with imaginary frequen- 
cies. This procedure makes it possible to first compute 



Wgg'(<1) w ) on ly f° r irreducible points of the first Bril- 
louin zone (BZ). We then determine easily the screened 
interaction for all k-points in the Brillouin zone using 
symmetry properties. 



3. Self Energy 

The self energy S is the key quantity of any GWA cal- 
culation. As previously noticed, we have chosen to avoid 
plasmon-pole models and compute S with the uj depen- 
dence of the screened interaction W within the RPA. 

First, we split the integral of the selfenergy into a bare 
exchange or Hartree-Fock contribution T, x and an en- 
ergy dependent contribution Y^ c (oj) which describes self- 
energy effects beyond T, x . The matrix elements of the 
self-energy are now given by the sum of 



(^i^i* kn ) = -|EEE 



\Mg n 



(k,q) 







m G 



|q + G| 



(12) 
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where the summation is over occupied states, and 
<* kn |£ c H|* k „) 

= ^EEE [Mg n (k, q)]* M£,"(k, q) 

q GG' rn 

x C G l G ,(k,q,c) (13) 

with 

Cg G/ (k,q jW ) 

2?r 7 lj + - e m (k - q) + i<5sgn(e m (k - q) - fi) 

where W c is defined as W c = W — v, with v being 
the bare Coulomb potential. To evaluate this integral 
directly on the real axis one should compute W c for 
many points u>' since the shape of W c along the real 
axis is rather ragged. Even though this has been done 
by some authors*^ we choose to avoid this difficulty 
by using the fact that W c is well behaved along the 
imaginary axis. In the present work, we have performed 
this integral using two different methods: 

In the first one, the contour of the frequency integral 
(|14f> is deformed in a way to obtain an integral along the 
imaginary axis plus contributions from the poles of the 
Green's function. In this case, we obtain the following 
expression: 



TABLE I: Calculated quasiparticle energies of silicon for some 
points in the Brillouin zone with our two different implemen- 
tations. The results are in good agreement with each other. 
In the last line, the minimum band gap E g is presented. 

First method" Second method 6 



T lv -11.85 -11.87 

r 25 < v o.oo o.oo 

Tigc 3.09 3.09 

r 2 / c 4.05 4.06 

Xi„ -7.74 -7.68 

X 4 „ -2.90 -2.91 

Xi c 1.01 1.03 

14X 4c 10.64 10.59 

L 2 , v -9.57 -9.50 

Li„ -6.97 -6.90 

L 3 /„ -1.16 -1.17 

Li c 2.05 2.03 

L 3c 3.83 3.83 

E g 0.92 0.90 



"Contour deformation method 
''Analytic continuation method 



analytically continued to the real axis by fitting it to the 
following Pade form 



C GG ,(k,q, W ) 
1 



u - e n (k - q) 



t/i ^"< G '(^-")(_ en(k _ q))2+ ^, 

±Wg G ,(q, ±(w - e„(k - q)))0(±(w - e„(k - q))) 
x 6(±{ui - M))0(±(e„(k - q) - »)) 

The first term represents the contribution along the imag- 
inary axis and is evaluated by Gaussian quadrature. 
The second is from the poles of the Green's function 
and its computation is done by fitting values of W 
at ±(w — e„(k — q)) from values on a given mesh of 
frequencies*^. Here /i denotes the Fermi level in the LDA 
and w" is defined to be real. This method is similar to the 
one used by Aryasetiawan for the implementation of the 
GWA based on the LMTO method in the atomic sphere 
approximation 7 (ASA), and within the GWA of Kotani 
and coworkers based on the full-potential linear muffin- 
tin orbital (FP-LMTO) methodi 7 . The reader can find 
more details about this integration procedure in Refs. 
7 17. Similar work has been also carried out by Bechst- 
edt and coworkers^* as well as by Fleszar and Hankc 21 
starting from a pseudopotential approach. 
In our second implementation, which is similar to that of 
Ref.0 we evaluate the matrix elements of the correla- 
tive part of the self-energy (\l/ m k|S c '(a;)|\E'ik) for a set of 
imaginary frequencies iw, the resulting quantity is then 



P(z) 



a + a±z + a 2 z 



- ajqZ 



N 



b\z + b2Z 2 



b M z 



(15) 



where a; and bi are complex parameters that are deter- 
mined during the fit along the imaginary axis. Values 
of 5 for N and of 6 for M provided us with an accu- 
rate and stable fit. The same kind of continuation has 
also been applied with success to compute the dynamical 
response function;22i2i so we are confident of its relia- 
bility. The main difference between the work presented 
here and that of Ref. lfl9l) is that our code starts from an 
all-electron basis, so we are not using fast Fourier trans- 
forms to switch between real and reciprocal spaces and 
between time and frequency domains. Our expression 
(|15f) is also different, but we believe that this is of minor 
importance^. 

In both cases, the integration over the first Brillouin 
zone is done by the special-point technique^ The num- 
ber of bands as well as the number of G vectors in i|13|) 
is increased until the QP energies are converged. Simi- 
larly, the number of frequency points lo' for which W c 
is computed is increased until C GG , (k, q, uj) is well con- 
verged. The two different implementations allow us to 
check carefully our results, and as can be seen in Table [I] 
for the case of silicon, the QP energies are insensitive to 
the method used to compute the self-energy. 
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4- Treatment of the Coulomb divergence 

The last point we wish to discuss is an additional dif- 
ficulty which occurs when evaluating the self-energy by 
a summation of the q points over the full BZ. We can- 
not apply the special-point technique directly since the 
integrands have a 1/q 2 singularity for q — > as can be 
seen for example by putting G = in the expression 
of the exchange term given by (|12f> . The difficulty can 
be removed by adding and subtracting a term which has 
the same singularities as the initial expression, and which 
can be evaluated numerically and analytically. As a con- 
sequence, the integrals over the BZ are rewritten 

£ c(q) = [ G (q) - A n<$\ + A E F (<*)> ( 16 ) 



where F(q) is an auxiliary periodic function that diverges 
like 1/q 2 as q vanishes. The term is regular and can be 
evaluated by the special point technique whereas the last 
sum is evaluated analytically. For the exchange term, it 
is not difficult to evaluate A in l|16|) , but it gets more com- 
plicated for the correlative part of the self-energy (|13fl . 
The purpose of the offsetted T-point methodic is to avoid 
the evaluation of the quantity A, but still to be able to 
deal with the divergence. The main idea is to find a new 
mesh of points such that 



$>(q)=£F(q') 



(17) 



where the T-point is included in the old mesh q but not in 
the new one q': the T-point is replaced by other points 
(different from T) to construct the q' grid in order to 
fulfill Eq. (TTJ. Eq. i|16|) is therefore rewritten as: 

£ G(q) = J2 [ G (q) - A ^(q)] +Aj2n<i'l (is) 



Then we show by inspection that the term 
E q [G(q) - A F(q)] is equal to £ q , [G(q') - A F(q')] 
with a controlled error, Eq. I|18fl transforms to 



^G(q) = ^TG(q') 



(19) 



because the two terms which contain the function F(q) 
cancel out since they are evaluated on the same q' grid. 
We have therefore avoided the evaluation of the compli- 
cated A term in Eq. l|ltj|) . 

The remaining points to be addressed are the choice of 
the function F(q) and the number of additional points 
for the new mesh used to solve ()17|l. In our case, we write 



f(q) 



^2 exp ( 

G 



q + G 



and choose to add 6 points to the original q mesh in 
order to get the new mesh. Equation l|17|) is then solved 



TABLE II: Lattice constants a (in atomic units) and energy 
cut offs E cut (in Rydberg) used for our PAW calculations. The 
lattice constants are from Ref. IjjH unless stated otherwise. 

a Ecut 



Si 

SiC 

AlAs 

InAs 

NaH 

KH 



10.26 
8.24 

10.67 

11.41 
9.28 a 

10.83" 



20 
25 
20 
20 
40 
40 



"Reference I 



to provide us with the coordinates of the new 6 points in 
the BZ.^Si The computational cost is further reduced by 
finding equivalent points among those 6 points, and we 
end up with well behaved and easily evaluated BZ sums. 



III. 



NUMERICAL RESULTS AND DISCUSSION 



In this section we present our theoretical quasiparti- 
cle energies for the six materials studied in this paper, 
together with the available theoretical and experimen- 
tal results. In Sec. IIII Al we report our results for four 
semiconductors (Si, SiC, AlAs, InAs) of zinc-blende type 
structure, while Sec. IIII Bl is devoted to studying the al- 
kali hydrides NaH and KH in the rocksalt phase. Table 
Irfl presents the experimental lattice parameters and the 
energy cut offs E cut used for the final converged calcula- 
tions. 



A. Results for Si, SiC, AlAs, and InAs 

Silicon is probably the most carefully studied semicon- 
ductor, and several GWA results are available. Using 
silicon as a prototype will allow us to test our method 
by making careful comparisons with previous GWA cal- 
culations. As mentioned earlier, our code presents two 
different ways for calculating the self-energy. We there- 
fore test their accuracy for silicon in Table [I] We find 
that the results of the two methods are almost identical, 
showing that they are equally reliable for computing the 
self-energy. In particular, our implementation with the 
extrapolation procedure makes it possible to represent 
the full-frequency dependence of the self-energy with a 
small additional computational cost. Fig. shows the 
real and imaginary parts of E along the real axis for sil- 
icon at the r point for a wide range of frequencies. The 
agreement with previous work is excellent^ 

A special feature of our work is the possibility of 
obtaining the imaginary part of the self-energy (see 
Fig. nj, a task virtually impossible within the PIP 



approximatio] 



.57 



The spectral function which can be 



obtained directly from the self-energy 
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FIG. 1: Re{* km |E(w)|* km > and Im(* km |E(w)|* km > shown 
for the first 8 bands for silicon at k = 0. The zero of energy 
is at the center of the band gap. 



<tf km |AH|tf 



|Im(* km |S(w)|* km }| 



e m (k) - Re(* km |E(w)|* km }] 2 + pm(tf km |£(a;)|* km >] 



is of major interest since it can be used for the in- 
terpretation of experimental photoemission and inverse- 
photoemission spectra. As an example, the spectral func- 
tion of silicon at the T point is shown in Fig. |3 The 
sharp peaks correspond to QP excitations, while the in- 
coherent part of the function, the spectral background, is 
much complicated and could correspond to plasmon type 
excitations. 

The QP calculations have been performed using 256 k 
points in the full BZ. The size of the dielectric matrix 
defined in Eq. © is 137 x 137 for silicon and SiC, 169 x 
169 for AlAs, 181 x 181 for InAs. 200 bands were used 
for the sum over conduction states in Eq. 10 and for the 
sum over m in Eq. 1|13[) . Due to the smoothness of the 
integrand along the imaginary axis, 11 points are found 
sufficient to obtain well converged quantities. An energy 
step of 1.5 eV is used for the part of Eq. i|15|) which 
corresponds to the poles of the Green's function. Using 



FIG. 2: Spectral function (* km |A(oj)|* km ) for the first 8 
bands for silicon at k = 0. The zero of energy is at the center 
of the band gap. The sharp peaks are QP poles, their weights 
correspond to the factor Z defined in Eq. 0] 



this energy step we determine an energy grid which we 
use to produce an accurate fit to W^ G ,(q, ±(cj — e„(k — 
q))) for the different points ±(w — e„(k — q)). All these 
high values of the parameters ensure the convergence of 
the QP energies to within 0.05 eV. 

Tabic IIIII shows the excellent agreement of our results 
with two other all-electron GWA implementations of the 
QP energies of silicon. From this table it seems that, at 
least for Si, the overall difference between the RPA and 
the PIP results is smallw Nevertheless, a discrepancy 
of as much as 0.18 eV for the energy of Li„ is obtained. 
It seems then, at least for Si, the PIP model overesti- 
mates only slightly the differences between the energy 
levels within the GWA. 

Table IIVI compares the calculated QP energies for 
3C-SiC (also known as /3-SiC), AlAs and InAs with 
experimental data as well as with pseudopotential-GWA 
(PP-GWA) calculations. The band gaps are given at T, 
X, and L and are underlined in this table. These studies 
2are motivated by the fact that SiC is a material of high 
current technological interest and that InAs is predicted 
to be metallic in the LDA, whereas AlAs is used as a 
test case. A general trend of our implementation is 
that the agreement with experiment as well as with 
PP-GWA results is not perfect, as also found by other 
implementations based on all-electron methodsji^i^ 
reporting differences up to 0.4 eV. In particular, the 
largest difference occurs for the T' 25v — > Xi c transition. 59 
We showed in a previous studyi 5 - for the case of silicon 
that these differences are mainly traced back to differ- 
ences between the exchange-correlation matrix elements 
obtained by the two methods. We believe that this can 
be extended to other materials as well, since it seems to 
be a general feature^ of all-electron GWA calculations. 

In fact, the difference between all-electron and PP 
based GW calculations is not surprising since the use 
of pseudo-wave functions for evaluating matrix elements 
of a general operator produce results that may not be 
sufficiently precise, because the pseudo-wave functions 
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TABLE III: Selected energy eigenvalues, in eV, at P, X and L for Si. Our results are compared with two other all-electron 
implementations of the GW method and with experimental results. Data in parentheses are results when the denominator of 
the Green's function is updated with QP energies. Data in the last line correspond to the minimum energy gap E g . 
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0.61 


0.65 


1.01 


(1.08) 


1.10 


1.14 


1.25 d 




10.02 




10.64 


(10.72) 


10.74 






T ' 


-9.63 


-9.63 


-9.57 


(-9.60) 


-9.66 


-9.92 


-9.3±0.4 


Liu 


-6.99 


-6.98 


-6.97 


(-7.00) 


-7.15 


-7.31 


-6.7±0.2 




-1.19 


-1.19 


-1.16 


(-1.17) 


-1.24 


-1.26 


-1.2±0.2 


Lie 


1.44 


1.43 


2.05 


(2.11) 


2.08 


2.15 


2.1 9 ,2.4±0.15 


L3c 


3.30 


3.35 


3.83 


(3.90) 


3.92 


4.08 


4.15±0.1 h 


E a 


0.55 


0.52 


0.92 


(0.95) 


0.97 


1.01 


1.17 



"Reference Y^a 

Reference [13 
c Unless noted, 
36 



^Reference 
E Referencc 



■^Reference 4£ 



9 Reference 
^Reference 



Rcf : 



.39 



41 
42 



are constructed to reproduce the all-electron wave 
functions only in the interstitial region, i.e., outside the 
atomic spheres. This is however a good representation 
for studying properties that depend only on the behavior 
of the wave function in the bonding region. An error is 
then introduced in any PP-GWA calculation. This error 
seems to fortuitously have a tendency to improve the 
agreement with experiment, explaining the exceptional 
success of the PP-GWA. 

For InAs, the incorrect metallic behavior obtained 
within the LDA is corrected by our GWA calculation. 
The GWA produces the true semiconducting state as 
given by experiment. Since we did not account for the 
spin-orbit coupling in our calculation, we have simply 
averaged out the spin-orbit-split experimental values to 
make the comparison with our work possible. 



B. Results for the alkali hydrides NaH and KH 

Alkali hydride materials exhibit a structural phase 
transition from the Bl (NaCl structure) to the B2 (CsCl 
structure) type structure under hydrostatic pressure. A 
number of studies have been performed to understand the 
equation of state of these materials, 36,37 ' 48 ' 49 as well as 
the possibility of an insulator-metal transition^SLSi how- 



ever only few studies have been published about the elec- 
tronic structures4242iSi In these materials, hydrogen be- 
haves as an H ~ ion, leading to partly ionic materials with 
a larger band gap than the studied 'sp' semiconductors. 
It is, therefore, of interest to know whether the GWA is 
capable of producing such large band gaps. In this study 
we are only concerned with the determination of their 
QP energies for the rocksalt crystallographic structure. 
In Fig. [3]we report the QP band structures of NaH and 
KH along some high-symmetry directions, 61 and present 
a detailed overview of their LDA and QP energies in Ta- 
ble In the following we detail and compare our re- 
sults with existing experimental and calculated results: 
For NaH, the authors of Ref. |48| reported an LDA cal- 
culation with an improved LMTO method in its atomic 
sphere approximation, however, from their band struc- 
ture plot we estimated that their band gap is only about 
2.7 eV and is direct at the L point. This disagrees with 
our calculation, since at the LDA level we found an indi- 
rect ba nd g ap of 3.39 eV from W to L. This could be due 
to Ref. [^Js use of a different value of the lattice param- 
eter of 8.90 a.u. A calculation by Kunz and Michish 49 
based on the electronic polaron model produced a direct 
band gap at X of 1.52 eV, a value too low for these ionic 
materials, so that it should be taken only at a qualitative 
levels However, our LDA results are in full agreement 
with results of Ref. l52t where an indirect band gap from 



8 



TABLE IV: Quasiparticle energies in eV at P, X and L for SiC, AlAs, and InAs. Data in parentheses are results when the 
denominator of the Green's function is updated with QP energies. Our results are compared with PP-GW method and with 
experimental results (minimum band gaps are underlined). 

LDA GW Expt" 





Present 




Present 


ppb 




SiC 












^25v — ► ri5 c 


6.25 


6.41 


7.32(7.45) 


7.35 






1.28 


1.31 


1.80(1.89) 


2.34 


2.39 


T*25v ~ * Lie 


5.34 


5.46 


6.45(6.56) 


6.53 


4.2 


AlAs 












^2&v ~~ * ri5c 


1.94 


1.77 


2.72(2.79) 


2.75 


3.11 c 




1.32 


1.20 


1.57(1.65) 


2.08 


2.24 


r25i> — * Lie 


2.06 


1.89 


2.73(2.80) 


2.79 


2.49 c ;2.54 


InAs 












T25t) — * Il5c 


-0.07 


-0.39 


0.46(0.49) 


0.59 


0.60 


r25i> ~~ * Xi c 


1.48 




1.57 (1.61) 


2.10 




T25J) ~~ * Lie 


1.05 




1.54 (1.58) 


1.52 


1.74 



"Unless noted Ref l35l For InAs, data have been averaged to 
account for the neglect of spin-orbit coupling in our case. 

6 For SiC, PP results are from Ref. [£| for AlAs from Ref. EJ 
for InAs from Ref. l4 r J and averaged to account for the neglect of 
spin-orbit splitting. 

c Reference |4|i 
Reference l47t 
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l r x w K 



FIG. 3: Calculated LDA (doted lines) and GW (full lines) 
electronic band structures of NaH and KH along some high- 
symmetry directions. These calculations are performed using 
the parameters reported in Table ITT1 

W to L of about 3.3 eV was found. Table IVl shows that 
the energy level differences of NaH are substantially in- 
creased by the use of the GWA compared with the LDA 
results. In particular, the minimum band gap is 5.68 eV 
within the GWA, whereas it's only 3.39 eV within the 
LDA. An other interesting point is that we found that 
the most used scissors-operator shift seems not to apply 
for the computation of optical properties of NaH. This 
is because the band gap is increased by 1.66 eV at the 
L point, and as much as 2.34 eV for the direct transi- 
tion at L point. A GWA calculation of the QP energies 
in the full BZ is then required for the study the optical 
properties of NaH. Table IVl shows also our LDA and QP 
results for potassium hydride. Our LDA results are in 
good agreement with the calculation in Ref. |52|- We find 
that in both LDA and GWA calculations, the band gap is 
indirect from W to L, and could easily switch to a direct 
gap with a small lattice parameter variation since the 
valence band (Is state of hydrogen) is very flat, i.e., the 
valences states at the X, L, W, and K have the same 
energies within 0.2 eV accuracy. The previous remark 
about the non-validity of the scissors-operator shift still 
holds here for the same reasons. The LDA band gap in- 
creases by an amount of 2.29 eV at L, and as much as 



TABLE V: Quasiparticle energies in eV at L, X, L, W, and K 
for NaH and KH. The last line shows the minimum band gap 
E g . Data in parentheses are results when the denominator 
of the Green's function is updated with QP energies. We are 
not aware of any experimental study concerning the electronic 
band structure of these materials. 





LDA 


NaH 

GW 


LDA 


KH 

GW 


ri„ 


-3.64 


-3 


39 


(-3.59) 


-2.07 


-2.02 


(-2.11) 




9.98 


11 


89 


(12.38) 


7.47 


9.81 


(10.43) 


1-25' c 


14.98 


16 


71 


(17.24) 


7.57 


10.16 


(10.74) 


Il5c 


15.72 


17 


87 


(18.11) 


9.33 


11.98 


(12.53) 


Xi v 


-0.13 


-0 


12 


(-0.12) 


-0.18 


-0.18 


(-0.20) 




4.00 


6 


23 


(6.66) 


4.91 


7.49 


(8.01) 




9.59 


11 


45 


(11.96) 


5.18 


7.92 


(8.47) 


-X5/ c 


9.95 


12 


15 


(12.60) 


8.76 


11.68 


(12.24) 


Xi c 


10.84 


12 


95 


(13.36) 


9.27 


12.31 


(12.86) 


Liu 


-1.03 


-1 


03 


(-1.06) 


-0.15 


-0.19 


(-0.20) 




3.39 


5 


68 


(6.11) 


3.18 


5.85 


(6.35) 


Lie 


10.24 


12 


13 


(12.59) 


6.86 


9.69 


(10.23) 




12.89 


15 


14 


(15.61) 


8.46 


10.49 


(11.18) 


L3 C 


15.42 


17 


34 


(17.77) 


10.98 


13.57 


(14.03) 


Wi, 


0.00 





00 


(0.00) 


0.00 


0.00 


(0.00) 


w 3c 


5.43 


7 


79 


(8.24) 


4.96 


7.60 


(8.13) 


w 2 , c 


7.75 


10 


08 


(10.55) 


6.70 


9.40 


(9.94) 


Wlc 


15.50 


17 


57 


(17.98) 


10.50 


13.03 


(13.55) 


w 3c 


17.08 


19 


16 


(19.57) 


10.73 


13.25 


(13.78) 


Kit, 


-0.13 


-0 


12 


(-0.12) 


-0.06 


-0.06 


(-0.07) 


K 3c 


4.71 


6 


96 


(7.41) 


4.84 


7.50 


(8.02) 


K ic 


5.84 


8 


11 


(8.56) 


4.94 


7.62 


(8.18) 


K 4c 


10.57 


12 


85 


(13.27) 


8.13 


10.72 


(11.38) 


K ic 


11.02 


13 


09 


(13.52) 


8.36 


10.94 


(11.54) 


E g 


3.39 


5 


68 


(6.11) 


3.18 


5.85 (6.35) 



2.71 eV at the L point. 

It is surprising to notice that, across the whole Bril- 
louin zone, the GW band gap shift of KH is larger than 
that of NaH despite that the LDA band gap of NaH is 
larger than that of KH. The reason for this puzzling large 
shift is that the screening of the Coulomb potential in KH 
is found to be less efficient than in NaH. Indeed, we have 
found that the RPA static dielectric function of NaH is 
3.43 much larger than that of KH which is about 2.62. 
The hybridization is also less strong for KH than for NaH, 
since the band width of hydrogen s-states is about 2.02 
eV for KH (much smaller than the 3.64 eV band width 
of NaH). As a consequence, the higher excited states are 
lower in energy, across the Brillouin zone, by about 2 to 
5 e V for KH than for NaH. We hope that our predictive 
results will stimulate experimentalists to perform photoe- 
mission or optical studies of these interesting materials. 
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IV. CONCLUSION 

We have presented an implementation of the GWA 
using the all-electron PAW method where the screened 
Coulomb interaction is obtained using the RPA dy- 
namical dielectric function. Thus we avoided the use 
of the plasmon-polc approximation. We have applied 
it to study the QP energies of Si, SiC, AlAs, and 
InAs and found that a precise comparison with other 
available theoretical and experimental results shows that 
sometimes the GWA can lead to noticeable discrepancies 
with experiments. Those discrepancies are generally not 
so pronounced in the pseudopotential GWA calculations 
using PIP models. We argued that the approximations 
used in the pseudopotential method have a tendency to 
fortuitously improve the agreement with experiment. 

We have presented detailed results for the first time 
for NaH and KH alkali hydrides, and showed that 
the GWA enhanced substantially the LDA band gaps, 
motivating further theoretical and experimental studies. 

Since our method can compute the imaginary part 
of the self-energy, we could then determine the QP 
lifetimes, a task not possible using the PIP approxima- 
tion. Further inspection of spectral properties as well 



as the computation of QP lifetimes will be presented 
in future work. The method is currently being applied 
to determine the excitation properties of LiH, and 
the results will be reported elsewhere^ Moreover, 
the use of symmetry and an efficient implementation 
make us confident that we will soon be able to study 
systems with a large number of atoms per unit cell, like 
surfaces or polymers. Furthermore, because we use a 
mixed basis-set in our implementation we could study 
systems with localized 'cf or '/' electrons with a reduced 
computational cost compared with methods based only 
a plane wave basis-set. 
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